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Abstract 

A systematic strategy for the calculation of density functionals (DFs) consists in coding informa- 
tions about the density and the energy into polynomials of the degrees of freedom of wave functions. 
DFs and Kohn-Sham potentials (KSPs) are then obtained by standard elimination procedures of 
such degrees of freedom between the polynomials. Numerical examples illustrate the formalism. 
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Existence theorems [1[ for DFs do not provide directly constructive algorithms. Fortu- 
nately, the Kohn-Sham (KS) method spares the construction of a "kinetic functional" 
and reduces energy and density calculations to the tuning of a local potential, VKs{f^)- Hence, 
a considerable amount of work has been dedicated to detailed estimates of electronic corre- 
lation energies and the corresponding KSPs, see for instance 0-[5[ . Many authors were also 
concerned with representability and stability questions, see for instance jcj and, for calcu- 
lations in subspaces, see and j^. For cases where the mapping between potential and 
density shows singularities, see ^ . For reviews of the rich multiplicity of derivations of DFs 
and KS solutions and their properties, we refer to [loj and and, for nuclear physics, to 

0- 

Local or quasi local approximations use the continuous infinity of values p(r), Vr, as the 
parameters of the problem. However, whether for atoms, molecules or nuclei, a finite number 
of parameters is enough to describe physical situations. For instance, Woods-Saxon nuclear 
profiles notoriously make good approximations, depending only on a handful of parameters, 
and it is easy to add a few parameters describing, for example, long tails and/or moderate 
oscillations of the density. High frequency oscillations are unlikely, for they might cost large 
excitation energies. The number of parameters to describe a nuclear density, therefore, can 
be restricted to maybe ~ 10 at most; situations with ~ 20 parameters are a luxury. For 
molecules, shapes are much more numerous, but a finite, while large number of parameters, 
truncating a list of multipoles for instance, still makes a reasonable frame. Practical DFs, 
therefore, can boil down to functions of a finite number of parameters. Integro-differential 
variations can then be replaced by simple derivatives. 

This Letter shows how informations about both the density and the energy can be recast 
into polynomials. This allows eliminations of part of the parameters. Further polynomial 
manipulations locate energy extrema. Only density parameters are left. The same method 
gives KSPs. Finally we offer a discussion and conclusion. 

A toy model of two particles illustrates the strategy. Consider the first 4 wave functions of 
the one-dimensional harmonic oscillator, ipo{r),. . . ,y93(r), and a Slater determinant $ made 
of positive and negative parity orbitals, 

= {U + w)ipo + Wip2, (j)- = {X + iy)Lpi + Z(p3, (1) 

where all 6 parameters, u, . . . ,z are real numbers. Parity ensures orthogonality. Normaliza- 
tions read, 

+ + = 1, x^ + y^ + z^ = l. (2) 
The density is a Gaussian, multiplied by a polynomial, 

p(r) = '!T^^e~''^{aer^ + a^r^ + a2r^ + ao), (3) 
with 3 independent coefficients, since the integral, 

dr p{r) = —ae + -a^ + -02 + ao = 2, (4) 
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equates to the particle number. Given a Hamiltonian if, the DF is defined by the constrained 
minimization jl3|, ^[p] = Min<j,=^p ($ \H\ $), where the constraint, $ ^ p, means that, 
given ag, 04, 02, (and hence ao,) the parameters, u, f , . . . , 2;, must ensure that, 

10+1^ + = 7r"^e-'''(a6r6 + a^r^ + a^r"^ + Oq). (5) 
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Insert Eqs. ([T]) into Eq. ([5]) and take advantage of the harmonic oscillator basis states. The 
density constraint, $ ^ p, then reduces into the 3 quadratic constraints, 
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Az = 3a6, Y ^ \ 3^^ ~ ^ 4 ' 
\Pluw — w"^ + + y"^ — \/Qxz + -z^ = (6) 

Recall there are 6 parameters, namely u,v, . . . ,z io handle such 3 constraints, Eqs. (EI), 
and also the 2 normalization constraints, Eqs. ([2]), and the energy. The latter, 77 = ($ \H\ $), 
is also a polynomial. Indeed, if if is a one-body operator, with given matrix elements, 
{ifi \H\(pj), the function 77 is quadratic with respect to the initial 6 parameters. If H contains 
two-body operators, then 77 has degree 4 with respect to the same 6 parameters, because of 
matrix elements, {^Pi^j \H \ ipk^i). 

For a model of maximum simplicity, consider the sum of 2 single-particle harmonic oscil- 
lators, without their zero-point energy, H = —(P/ {2drf) — (P / {2dr'^ + {r\ + r|)/2 — 1. Then 
the energy reads, 

■q = 2w'^ + + y"^ + ?>z^ . (7) 
A trivial elimination of 5 among the parameters, by means of the 5 constraints, gives, 

P(?7, = 33 - 12a2 + 4a^ - 1804 + 12a2a4 + 904- 

Gae + 360206 -l- bA^a^a^ + 99ag — IQv'^ — 24a6f ^ — 

ri{2Q + 4a2 + 604 + 42a6 - IQv"^) + ^rf = 0. (8) 

There remains to minimize rj with respect to that last degree of freedom, v, left after the 
elimination. Let a "precursor" be a polynomial resulting from the definition of the energy, 
Eq. ([7j), after all constraints have eliminated a maximum number of degrees of freedom. In 
this toy model, a precursor is the left-hand side, P(?7, v), of Eq. ([8]). Extrema of r] correspond 
to the cancellation of, 

dP/dv = 16v{2r]-3aQ-2). (9) 

That solution which directly relates rj to Og only is clearly spurious. The correct solution is, 
therefore, v = 0. Accordingly, the rigorous link between 77 and the parameters of the density 
is given by the polynomial equation, 

= P{ri, 0) = 33 - 12a2 + Aaj - ISa^ + 120204 + 9a^ - Ga^ 

+ 3602^6 + 540406 + 99o6 - (26 + 4o2 + 604 + 42o6)?7 + 9?7^ (10) 

This is not an open formula of the form, rj = F[o2, 04,06], but it provides roots for rj at 
any realistic degree of numerical accuracy. Incidentally, it gives excited energies as well as 



minima, a well known property [15| of DFs. 

For numerical tests we show in Fig. [T]a plot of a density, pi(r) = (0.37861 + 3.64114r^ — 
3.30773r'^ + 1.21686r^)e~'" /a/vt. Fig. [2] shows the r] roots of Eq. ([8]) when v takes on all 
values that induce non-complex roots. We took great care to verify that the same values of 
Tj are obtained if v is run as a parameter of Eqs. (jS]) and Eqs. ([2]) to obtain u, w, x, y, and 
z by a brute force numerical algorithm without the elimination process. 

Figs. E] and S show another density, p2{r) = (1.11314 - 0.0853536r2 + 0.10563r^ + 
0.453499r^)e"'' /y/n, and the r] roots obtained therefrom. In both cases of the density. 
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FIG. 1: The first density, pi{r), as described in the text. 
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FIG. 2: Energies rj compatible with the density of Fig. [H in terms of the free parameter v. 



and in every other case that we tested at random, the minimum of for f = is confirmed. 
It may be noticed that, because of a better centering and milder oscillations for the density 
seen in Fig. [3|, its pattern of energies, see Fig. |U lies lower than that seen in Fig. [2l 

A side issue arises here, and will also arise in all future models using this polynomial 
method, because a final minimization of rj must be performed within a convex domain of 
densities: what conditions the coefficients aj (and their resulting oq) must satisfy to maintain 
p(r) positive? This question was recently [l6| solved by means of the Sturm criterion, for 
positive functions having positive Fourier transforms. The criterion gives the number of real 
roots of a polynomial, and can be used to ensure that a polynomial has no real roots. In the 
present model, the minimum energy among all densities corresponds to that determinant 
made of ipo and with density, p(r) = (1 + 2r^)e~^ / V^- If we set, therefore, 02 = 2 and 
CL4 = Ctg = in Eq. (ITOl) . it is easy to verify that the root, r] = 1, is found, as expected from 
Eq. ([7]) with u = X = 1 and v = w = y = z = 0. 

Returning now to the "only kinetic" KS functional, we keep Eqs. ([2]) and (Q, but replace 
Eq. dZD by, 

u"^ + v'^ - 2y/2 uw + + + 3y^ - 2y/6xz + Tz"^ 



T 
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FIG. 3: The second density, /02(r), as described in the text. 
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FIG. 4: As for Fig. [21 but for the density P2{r) in Fig. [3l 



Elementary eliminations provide a "precursor" , 

Q{t, v) = 12352 - 480002 + 912a^ - 10560a4 + 2784a2a4+ 

2448a^ - 2136006 + 68880306 + 111600406 + 13329o^ - 8960r 
+ 128o2r + 1920o4r + 2208o6r + (-6144 + 1024o2+ 

3072O4 + 537606 + 4096r)t;2 + 2304r2. (12) 

The final equation for our "algebraic" brand of KS kinetic functional reads, 

= g(r, 0) = 12352 - 48OO02 + 912o^ - IO56O04+ 
27840204 + 2448o^ - 21360o6 + 688802O6 + III6O04O6+ 
13329o^ + (-8960 + I2802 + I92O04 + 2208o6)r + 2304rl 

The relation, t = [i] + 1) — u, that connects the kinetic energy to the potential one. 



(13) 
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r^dr p{r) = (8 + 4o2 + I204 + 45o6) /16, 



(14) 
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and to the full energy, transforms Eq. (fT3ll into Eq. (fTOl) . 

Among generalizations of some interest, consider the radial density functional (RDF) 
theory [l7|, valid for nuclei and/or atoms, isolated, described by rotation invariant Hamil- 
tonians. The constrained density minimization of energy isl, l3 returns isotropic densities. 



with radial profiles, p(r), ^ r < oo. Choose a positive measure, /i(r), having some physical 
relevance, and find a basis {pen} of orthonornal polynomials, to define a set of "raw orbitals" , 
<^n£m(r) = Yimi^r) Kr)Pen{r). Physical orbitals, 0^ = J^nim'^nem^nem, can be expanded on 
this basis, truncated at some maximum polynomial order. Furthermore, whatever the pro- 
ton, neutron or electron numbers, states $q made of such orbitals can be mixed to generate 
isotropic densities p(r). Any such p is a polynomial of r multiplying [p(r)]^. Simultaneously, 
the density, and the energy, will be polynomials with respect to the coefficients c"^^ and/or 
the coefficients of configuration mixing, \E' = Ylg^q^q- There will therefore remain free 
parameters once a maximum number of parameters has been eliminated between the energy 
polynomial and the polynomials describing the density fit constraints. These free parameters 
stay in what we called a precursor polynomial. Energy minimization by manipulations of 
such a precursor requires again standard algebraic operations only. In a forthcoming paper, 
we shall describe a model more realistic than that used for this Letter. 

This strategy gives directly the full DF, hence making the KS detour useless in principle. 
If spurious solutions 0,[8| pop up, an analysis for their detection is easy, because the polyno- 
mial equation for the energy can only create a finite number of candidate solution branches. 
The KS approach may remain useful to verify n- and w-represent ability. A constructive 
derivation of KSPs is available, indeed. For instance, truncate again the harmonic oscillator 
basis at a maximum polynomial order n for single particle wave functions and let V be the 
projector upon the resulting, finite dimensional subspace for a system of fermions, with 
their Hamiltonian H, or rather now, VHV. Any density in the subspace shows, besides 

2 

its exponential factor, e'^"' , a polynomial of maximum order 2n. Given the kinetic energy 
operator T, choose a local potential fo(r), hence a one-body operator Vq = Yld=i'^o{^i)^ 
hence again a one-body Hamiltonian Hq = T + Vq, so that the ground state of VHqV, a 
Slater determinant $0; be non degenerate and providing an approximate density po for the 
system. For any density p in the subspace, the difference, Ap = p — po, with J Ap = 0, 
remains the product of e'"^^ by a polynomial of maximum order 2n. (Here and in the 
following, the integral sign, J, means J r'^"^ dr depending on the d-dimensional problem un- 
der consideration.) Expand Ap in a basis of orthonormal polynomials >S'/3(r), "constrained 
by vanishing averages" [l8l-20|. Ap(r) = ^/3(?"), where ^/3(r) = e~'^^^ Sp{r). Again, 

given a determinant $ with the parameters c"^^ of its orbitals, or given a correlated state, 
\1/ = ^gC*g$q, the constraints, $ ^ 6/3 or \E' ^ 6/3, are polynomials of the parameters. 
Given Hq, the polynomial method provides a polynomial /C(fi;; 61, . . . , 62^) ioT a reference 
functional, such that the lowest root of the equation, /C = 0, represents the constrained 
minimum, k' = Minij,^hi,...,b2„{^\Ho\^) , for the determinants in the subspace. In the same 
way, given the full H, the method provides a polynomial £{ri;bi, . . . , b2n), the lowest r] root 
of which is the constrained minimum, t]' = Mm^^bi,...,b2„{'^\H\'^) , for correlated states in 
the subspace. Then it is trivial to derive from JC and S a polynomial, Q{u; 61, ... , b2n), for 
the difference, u = t] — n. The diagonalization of VHV then reads. 

Ok' dco' ^ , . 

With the ratio, = —{dfl/db0)/{dfl/du), representing du' /dbp, define the one-body, local 
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potential, v^ir) = X]^=i "^z? ^/^l'^)- Let $ be the ground state of V Hq + J^iLi'^Ai^i) 
Notice that {^iV^fsV]^) = ($|,^^|$). Then the energy E of $ has derivatives, 

dEldvp = J {Ap + po) ^0 = bp + bpo, (16) 

because of the orthornomahty of the S'^'s under the weight, e"^'^^^ The numbers, 6^0 = 
/ pQe~"'^ Sjs, are easily pretabulated. The quantities, and (6/3 + 6^0)) are Legendre con- 
jugates, and, moreover, d/d{hp + 6^0) = d/dhj3. The conditions, Eqs. ( |T5i) . read as the 
diagonalization for a determinant $ with the same density p as that of the eigenstate \E' of 
VHV. The potential, V{vq + v^)V, is a KSP valid for the subspace. 

Rather than coefficients, we also fitted local values, p(rj), i = 1, . . . ,n, for even profiles. 
Such fits are also polynomial constraints with respect to the parameters. An optimization 
of the choice of such points makes a side issue, of some interest. But, at present, we feel 
that the simplest algebra results from fitting density polynomial coefficients. This is a highly 
nonlocal parametrization of p that deviates from the quasi-local tradition of the field. Other 
parametrizations, such as by moments of p from a Fourier transform, deserve consideration. 
In every case, our unconventional parametrization of p, in the frame of a polynomial algebra, 
creates a new zoology of DFs. Nothing of this zoology is known to us, but its interest is 
obvious, since manipulations of polynomials and properties of their roots, including bounds, 
are basic subjects. The number of available, exactly solvable models is huge, both at zero 
and finite temperature. It is limited only by computational power. For nuclei or atoms. 
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the models will be "radial" [17[, hence somewhat simple, although density tails adapted to 
Coulomb potentials, even screened ones, might demand many components. For electrons 
in molecules or extended systems (metals, thin layers, etc.), however, a necessary algebra 
of polynomials of 2 or 3 variables will burden the models. Anyhow, one can always test 
whether functionals from "smaller" models may remain valid for "larger" ones. Asymptotic 
properties might guide towards more traditional derivations of DFs. The models allow 
comparisons between the KS and the true kinetic energies of correlated systems. They also 
provide DF terms for those correlation energies due to interactions. 
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